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^ : ABSTRACT 
> ■ 

^ : 

■ We use three-dimensional (3D) numerical magnetohydrodynamic simulations to fol- 
QQ . low the evolution of cold, turbulent, gaseous systems with parameters chosen to repre- 

, sent conditions in giant molecular clouds (GMCs). We present results of three model 

I cloud simulations in which the mean magnetic field strength is varied {Bq = 1.4 — 14 

"Th I for GMC parameters), but an identical initial turbulent velocity field is introduced. We 

Qh' describe the energy evolution, showing that (i) turbulence decays rapidly, with the tur- 

O ' bulent energy reduced by a factor two after 0.4-0.8 flow crossing times (~ 2 — 4 Myr for 

' GMC parameters), and (ii) the magnetically supercritical cloud models gravitationally 

■ collapse after time ^ 6 Myr, while the magnetically subcritical cloud does not collapse. 
J> . We compare density, velocity, and magnetic field structure in three sets of model 

^ , "snapshots" with matched values of the Mach number ~ 9, 7, 5. We show that the 

H ' distributions of volume density and column density are both approximately log-normal, 

■ " " ' with mean mass-weighted volume density a factor 3 — 6 times the unperturbed value, 

but mean mass- weighted column density only a factor 1.1 — 1.4 times the unperturbed 
value. We introduce a spatial binning algorithm to investigate the dependence of kinetic 
quantities on spatial scale for regions of column density contrast (ROCs) on the plane 
of the sky. We show that the average velocity dispersion for the distribution of ROCs is 
only weakly correlated with scale, similar to mean size-linewidth distributions for clumps 
within GMCs. We find that ROCs are often superpositions of spatially unconnected 
regions that cannot easily be separated using velocity information; we argue that the 
same difficulty may affect observed GMC clumps. We suggest that it may be possible 
to deduce the mean 3D size-linewidth relation using the lower envelope of the 2D size- 
linewidth distribution. We analyze magnetic field structure, and show that in the high 
density regime n/fj ^ lO^cm^^, total magnetic field strengths increase with density 
with logarithmic slope ~ 1/3 — 2/3. We find that mean line-of-sight magnetic field 
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strengths may vary widely across a projected cloud, and are not positively correlated 
with column density. We compute simulated interstellar polarization maps at varying 
observer orientations, and determine that the Chandrasckhar-Fcrmi formula multiplied 
by a factor ~ 0.5 yields a good estimate of the plane-of sky magnetic field strength, 
provided the dispersion in polarization angles is ^ 25°. 



1. Introduction 

Since the identification of cold interstellar clouds in radio molecular lines, observational cam- 
paigns in many wavelengths have provided an increasingly detailed and sophisticated characteriza- 
tion of their structural properties. These clouds are self-gravitating entities permeated by magnetic 
fields and strongly supersonic turbulence; the observational properties of giant molecular clouds 
(GMCs) are summarized, for example, by Blitz (1993); Williams, Blitz, & McKcc (2000); Evans 
(1999). Although it has long been appreciated by theorists that turbulence and magnetic fields 
must play a decisive role in cloud dynamics (e.g. Mestel & Spitzcr (1956); Shu et al. (1987); McKee 
et al. (1993); Shu et al (1999); McKee (1999)), much of the theoretical emphasis has been on 
evolutionary models in which the effects of turbulent magnetohydrodynamics (MHD) is modeled 
rather than treated in an explicit fashion. 

Recent advances in computer hardware and development of robust computational MHD al- 
gorithms have now made it possible to evolve simplified representations of molecular clouds using 
direct numerical simulations. Fully nonlinear, time-dependent, MHD integrations can test theoret- 
ical ideas about the roles of turbulence and magnetic fields in cloud evolution, and also make it 
possible to investigate how turbulence affects the structural properties of clouds. Progress in the 
rapidly-developing field of simulations of GMC turbulence is reviewed by, e.g., Vazquez-Semadeni 
et al (2000). 

This is the fourth in a series of papers (Gammie & Ostriker 1996; Stone, Ostriker, & Gam- 
mie 1998; Ostriker, Gammie, &; Stone 1999) [Papers I-HI] investigating the dynamics of turbulent, 
magnetized, cold clouds using direct numerical simulations. The previous papers presented several 
results on energetics and overall cloud evolution. They showed that (1) MHD turbulence can delay 
gravitational collapse along the mean magnetic field in one-dimensional models since dissipation is 
slow (Paper I); however (2) in higher dimensional models dissipation occurs on the flow crossing 
timescale tj (Paper II); as a consequence, (3) the fate of a cloud depends on whether its mass-to- 
magnetic flux ratio is subcritical or supercritical, independent of the initial turbulent excitation, 
provided that turbulence is not steadily driven (Paper III). 

Some important astrophysical implications of these results are that (1) star formation in tur- 
bulent clouds may be initiated rapidly, essentially on a flow crossing timescale; and (2) models that 
rely on slowly dissipating turbulence to support GMCs against collapse do not appear to be viable. 
One is still faced with the problem of avoiding the excessive Galactic star formation rate that 
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would result from the collapse and fragmentation of the whole cold component of the ISM within 
its gravitational free-fall time (comparable to its flow crossing time; Zuckerman & Palmer (1974)). 
This requires either (a) limitation of the star formation rate in individual clouds (if self-gravitating 
clouds are long-lived after formation), or (b) limitation of the lifetimes of self-gravitating clouds. 
Both of these effects may be important. Processes that contribute to limiting star formation rates 
in individual clouds include turbulent feedback from star formation, transmission of turbulence 
from the larger-scale ISM, or a large (subcritical) mean magnetization of clouds. The first two of 
these processes, together with destabilizing environmental factors such as enhanced galactic shear 
outside spiral arms, contribute to limiting lifetimes of individual clouds. 

Other workers have independently used simulations to deduce the same results about the 
rapidity of turbulent dissipation under likely GMC conditions (Mac Low et al 1998; Mac Low 
1999; Padoan &; Nordlund 1999). Similar conclusions have also been reached concerning ongoing 
turbulent driving and the potential for star formation to be initiated on a rapid timescale (see also 
Ballesteros-Paredes, Hartmann, & Vazquez-Semadeni (1999); Elmegreen (2000); Klessen, Heitsch, 
& Mac Low (2000)). 

In addition to studying cloud evolution, our previous work also investigated structural prop- 
erties of our model clouds. We found that density contrasts produced by turbulent stresses are 
compatible with the typical clump/interclump ratio estimated in GMCs (Papers I-III). We also 
found that velocity and magnetic field power spectra evolve to be comparable to power-law forms 
of Burgers and Kolmogorov turbulence, regardless of the driving scale (Paper I; see also Stone, 
Gammie, & Ostriker (2000)). Other workers have also studied the basic structural properties of the 
turbulent gas in compressible hydrodynamic and MHD simulations, concentrating on distribution 
functions of density and velocity (Vazquez-Semadeni 1994; Passot & Vazquez-Semadeni 1998; Scalo 
et al 1998; Nordlund & Padoan 1999; Klessen 2000), the ability of stresses to produce transient 
structure (Ballcstcros-Paredes, Vazquez-Semadeni, &: Scalo 1999), and power spectra and related 
functions (Passot, Vazquez-Semadeni, & Pouquet 1995; Vazquez-Semadeni, Ballesteros-Paredes, Sz 
Rodriguez 1997; Elmegreen 1997, 1999; Klessen, Heitsch, & Mac Low 2000; Mac Low & Ossenkopf 
2000). 

In this paper, we analyze decaying turbulence in self-gravitating cloud models with varying 
mean magnetization (i.e. mass-to-magnetic flux ratio). ^ We begin by briefly describing the 
energy evolutions of our models, which serve to conflrm our earlier results on turbulent dissipation 

times and the gravitational collapse criterion for magnetized clouds. We then turn to detailed 
structural investigations. We analyze the density, velocity, and magnetic field distributions in our 
models at those stages of evolution when the turbulent Mach number is comparable to that in 
large (~ 5 — 10 pc scale) clouds. Our goals are (1) to provide a basic description of structural 
characteristics and how they depend on input parameters; (2) to make connections between cloud 



^AU models reported here, and most models studied by other workers, impose the somewhat artificial constraint 
that the initial mass-to-flux ratio is spatially uniform. 
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models seen in projection and their true three-dimensional (3D) structure, so as to help interpret 
observational maps; and (3) to assess whether certain statistical properties of clouds can be used 
to estimate the mean magnetization. 

Wc divide our analysis of structure into three main sections. The first (§4) is a discussion of 
density structure. Numerical hydrodynamic and MHD simulations of supersonic, turbulent flows 
have shown that magnetic pressure and ram pressure fluctuations produce structures with large 
density contrast that appear to resemble analogous "clumpy and filamentary" structures in real 
clouds (e.g. Passot, Vazquez-Semadeni, & Pouquet (1995); Padoan k, Nordlund (1999); Ostriker, 
Gammie, & Stone (1999); Klessen (2000); Balsara et al (1999)). Studies of density maxima and 
their immediate surroundings ("clumps") show that many are transient, as indicated by comparable 
values for the kinetic energy and kinetic surface terms in the virial theorem (Ballesteros-Paredes, 
Vazquez-Semadeni, k, Scalo (1999), cf. McKee k, Zwcibcl (1992)). Clump properties in our turbu- 
lent cloud models will be examined in a companion paper (Gammie, Stone, & Ostriker 2000). 

Analysis of the correlations of overdensity - via clump studies or multi-point statistics - will 
be needed to characterize fully how a spectrum of self-gravitating condensations is established. 
This process is of great interest because it may ultimately determine the stellar IMF. A first step 
in understanding the effect of turbulence on density structure, however, is to examine one-point 

statistics. Here, we compute and compare the distributions of density and and column density in 
different cloud models. We consider both because volume densities can be inferred only indirectly 
from observations, whereas column density distributions can be obtained directly from surveys of 
stellar extinction to background stars (Lada et al (1994); Alves et al (1998); Lada, Alves, k Lada 
(1999)). 

The second structural analysis section (§5) considers the linewidth-size relation. Observations 
give differing results for the slope of this relation depending on whether the structures involved 
are clearly spatially separated from the surroundings (e.g. by a large density contrasts), or are 
identified as coherent regions in position-velocity maps. The former case yields relatively steep 
power spectra (Larson 1981; Solomon et al 1987); the latter case yields shallower power spectra 
and larger scatter (Bertoldi k McKee 1992; Wilhams, de Geus, k Blitz 1994; Stutzki k Glisten 
1990), and has led to the concept of moderate-density "pressure-confined clumps" within GMCs. 
We believe the different slopes are a consequence of different definitions of "clump." We use a simple 
binning algorithm to explore the scaling of kinetic properties of apparent clumps within projected 
clouds, and in particular to understand the consequences of projection effects for lincwidth-size 
relations for 2D areas and 3D volumes. We argue, consistent with the suggestions of some other 
workers (e.g. Adler k Roberts (1992); Pichardo et al (2000)) that it may be difficult to identify 
spatially-coherent condensations from observed position-velocity maps. 

The third structural analysis section (§6) considers the magnetic field. A topic of much interest 
in turbulence modeling is understanding how the magnetic field affects both the intrinsic d3mamics 
and the observable properties of a cloud. As shown in §3, a major dynamical effect of the magnetic 
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field is to prevent gravitational collapse in subcritical clouds. Because magnetic field strengths are 
difficult to measure directly, however, it is highly desirable to determine if more-readily observable 
structural properties of clouds could act as proxies for the magnetic field strength. With simulations, 
it is possible to make comparisons of different models in which the mean field strength is varied, 
but other key properties (such as the turbulent Mach number and power spectrum) are controlled. 

An important theme in our analysis in §4 and §5 is to test how the quantitative measures of 
structure depend on the mean magnetization, to evaluate the potential use of such measures as 
indirect magnetic field diagnostics. In §6, we analyze how more direct magnetic field diagnostics 
may be affected by cloud turbulence. We evaluate the distribution of total magnetic field strength 
as a function of density in cloud models with different mean magnetization. We also compute the 
distribution of mean line-of-sight integrated magnetic field (one-point statistic), which is relevant 
for interpreting Zeeman effect measurements of magnetic field strength. Finally, we study the dis- 
tribution of polarization directions in simulated maps of polarized extinction produced by turbulent 
clouds (one-point statistic). One of the earliest estimates of magnetic field strength in the interstel- 
lar medium (Chandrasekhar & Fermi (1953)) was based on the dispersion in polarization direction, 
using a simple one-wave description of the magnetic field. Wc update the Cliandrasckhar-Fcrmi 
(CF) estimate using our simulations as presumably more realistic descriptions of the magnetic field 
geometry. 

The plan of this paper is as follows: We start (§2) by describing our numerical method and 
model parameters. We then (§3) describe our results on energy evolution, confirming the previous 
results from 3D non-self-gravitating models on dissipation rates, and 2.5D self-gravitating models 
on the criterion for collapse. We present our structural analyses in §§4-6, and conclude in §7 with 
a summary and discussion of these investigations. 

2. Numerical method and model parameters 

We create model clouds by integrating the compressible, ideal MHD equations using the ZEUS 

code (Stone & Norman 1992a, b). ZEUS is an operator-split, finite-difference algorithm on a stag- 
gered mesh that uses an artificial viscosity to capture shocks. ZEUS uses "constrained transport" 
to guarantee that V • B = to machine precision, and the "method of characteristics" to update 
the magnetic field in a way that ensures accurate propagation of Alfvenic disturbances (Evans & 
Hawley 1988; Hawley & Stone 1995). The solutions are obtained in a cubic box of side L with 
grids of 256^ zones, which permits spatial resolution over a large dynamic range at manageable 
computational cost. We apply periodic boundary conditions in all models. The simulations were 
run on an SGI Origin 2000 at NCSA. 

For the energy equation, we adopt an isothermal equation of state with sound speed Cg. In 
the absence of a fully time-dependent radiative transfer, this represents a good first approximation 
for the gas at densities higher than the mean - comprising most of the matter - for conditions 
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appropriate to molecular clouds (see discussion in Paper III, and also Scalo et al (1998)). 

The gravitational potential is computed from the density using standard Fourier transform 
methods. The A; = components of the density are not included in the solution due to the 
periodic boundary conditions. Rather than the usual Poisson equation, the gravitational potential 
therefore obeys V^^g = 47rG(p — p), where p = M/L^ is the mean density (mass/volume in the 
box). 

The initial conditions are as follows: We start with uniform density, a uniform magnetic field 
Bq = Bqx, and a random velocity field 5v. As in our earlier decay models (Papers I-III), 5v is a 
Gaussian random perturbation field with a power spectrum \5vk\^ oc subject to the constraint 
V • (5v = so that the initial velocity field is non-compressive. This power spectrum is slightly 
steeper than the Kolmogorov spectrum (|(5vfep oc k~^^/^) and matches the amplitude scaling of the 
Burgers spectrum associated with an ensemble of shocks (but differs from Burgers turbulence in 
that the initial phases are uncorrelated) . 

In configuration space, the velocity dispersion of the initial conditions averaged over a volume 
of linear size R increases as oc B}/"^ . This spectrum is comparable to the spectrum inferred for 
large-scale cold interstellar clouds (e.g. Larson (1981); Solomon et al (1987); Hcyer & Schlocrb 
(1997)) and the spectrum that naturally arises from the evolution of compressible turbulence that 
is either decaying or is driven over a limited range of scales (Stone, Gammie, & Ostriker 2000). We 
use an identical realization of the initial velocity field for all of the models, so that initial states of 
the simulations differ only in the strength of the (uniform) mean magnetic field. 

This paper considers three different simulated cloud models. All are initiated with kinetic 
energy = lOOpL^c^, corresponding to initial Mach number M. = o"„/cs = 10\/2. For the purposes 
of comparison with observations, we shall use a fiducial mean matter density (i.e. corresponding 
to the total mass divided by total volume) nu^ = lOOcm^'^ and isothermal temperature T = lOK 
in normalizing the local simulation variables of our models to dimensional values. The velocity 
dispersion in physical units is given by = 0.19 X M kms-^(r/10K)V2^ go that the initial value 
is ay = 2.7 kms-i(r/10K)V2. 

The models differ in their initial magnetic field strength, parameterized by /3 = c^/v^q = 
Cg / {Bq / Airp) , with physical value given by 

B„ = 1.4xr'/VG(^)"V ""^^ V'\ (1) 
^ ^ \WK J VlOOcm-3/ ^ 

We run a "strong field" model with /5 = 0.01, a "moderate field" model with (3 = 0.1, and a "weak 
field" model with /3 = 1. For characteristic fiducial densities and temperatures of molecular clouds 
(T ~ lOK, ~ 100cm~3), the corresponding uniform magnetic field strengths are 14, 4.4, and 
1.4 /xG. Of course, the evolved fields are spatially nonuniform and can differ greatly from these 
initial values (see §6), although the mean magnetic field (i.e. the volume-averaged value or A; = 
Fourier component) is a constant Bqx in time. The values oi (3 - half the ratio of the gas pressure 
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to the mean field magnetic pressure - are proportional to the square of the mass-to-magnetic-fiux 
ratio in the simulation box; this ratio cannot change in time. 

We may identify several different, physically significant timescales in the model evolution. The 

sound crossing time, tg = L/cg, is fixed due to the isothermal equation of state. Another important 
timescale is the flow crossing time over the box scale L, tf = = 9.8Myr x (^i^^ ( kms-i ) ^ " 
Because the turbulence decays (i.e. Ai decreases), the instantaneous flow crossing time increases 
relative to the sound crossing time as tf = tg/M. Where wc relate tf and i^, we use the Mach 
number associated with the imtm? turbulent velocity dispersion, (T^/cg = 14.1, such that tf = Q.Qltg. 

This paper concentrates on structures that form as a consequence of turbulence, before self- 
gravity becomes important. However, we also use the present models to test our previous results 
from lower-dimensional simulations (Paper III) on the differences in the gravitational collapse times 
with strong and weak mean magnetic fields Bq. It is therefore useful to define a gravitational 
contraction timescale 

In the absence of self-gravity, the unit of length L defining the linear scale of the simulation 
cube would be arbitrary. In a self-gravitating simulation, an additional parameter must be chosen 
to represent the relative importance of gravity and thermal pressure forces to the evolution. A 
useful dimensionless measure of this is tg/ts', in all the models considered here this ratio is 1/3. A 
more transparent way of stating this is that there are three thermal Jeans lengths Lj = Cgij^ IGpf-l"^ 
across a box scale L? 

The three simulations described herein differ in the relative importance of magnetic and grav- 
itational forces to their ultimate evolution. As described in Paper III, a cloud with constant 
mass-to- flux ratio is super- or sub- critical if tg is smaller or larger than ttL/va, respectively. A 
supercritical (subcritical) cloud has a ratio of mass-to-magnetic-flux greater (smaller) than the crit- 
ical value, l/(27rG^/^). Subcritical clouds can collapse along the field but not perpendicular to the 
field ( "pancake" ) ; in the nonlinear outcome the peak density would be limited by the thermal pres- 
sure. Supercritical clouds can collapse both parallel and perpendicular to the field, with unlimited 
asymptotic density. The three models discussed here have tgVA_/{7rL) = 0.11,0.34, and 1.1. Thus, 
the strong-field model is subcritical and the other two models are supercritical. The results on 
long-term gravitational evolution reported in §3 confirm the expected differences between super- 
and sub- critical clouds under the condition that turbulence secularly decays. 

Since self-gravity is weak for the first portion of the evolution in our simulations, the freedom 
of normalization of L that applies to non- self-gravitating models also effectively applies during 
this temporal epoch. In particular, the structural analyses of §§4-6 are performed at stages of the 



^Of course, the presence of strong turbulence makes the classical, linear Jeans stability analysis inapplicable; the 
velocity field is in the nonlinear regime from the first instant. 
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simulations' evolutions for which the kinetic energy is at least five times as large as the components 
of the gravitational energy Eq associated with the fluctuating density distribution. Because of 
our periodic boundary conditions, the gravitational energy associated with the mean density (i.e. 
the k = Fourier component) is not included in Eg- In order of magnitude, the value of this 
lowest-order gravity is ~ GM'^/L = Mc^Tr{L/ Lj)^ , which equals 28MCg for our present L/Lj = 3 
models. This energy is in the middle of the range of kinetic energies for the = 9, 7, 5 snapshots 
we analyze. Thus, as for observed clouds (e.g. Larson (1981); Myers & Goodman (1988)), the 
lowest-order gravitational energy in these snapshots would be comparable to the kinetic energy. 

A useful reference length scale may be obtained by combining the well-known observational 
relations between velocity dispersion, mass, and size for GMCs (see Paper III). The characteristic 
outer linear size scale Lobs for observed clouds scales with Mach number M according to 



Because the observed scale is proportional to the Mach number, the flow crossing time for observed 
clouds is independent of M , and given by 



For observed clouds, the flow crossing time and gravitational contraction time are proportional, 
with tf^obs ~ 0-59tg. 

Since the turbulence (and therefore M) decays in our models, they are comparable in their 
kinetic properties to increasingly small clouds as time progresses. For example, using the relation 
(3), the observational scale associated with the initial models with A4 = 14.1 would be i^o6s,imt = 
16pc X (r/10K)^/2 (nH2/100cm-3)"^/^ At this size scale, the corresponding sound crossing time 



properties of model snapshots in which ~ 9, 7, and 5; observed clouds of linear size scale ~ 10, 
8, and 6 pc, respectively, have kinetic energies corresponding to those of the model snapshots. To 
the extent that gravity may be unimportant for much of the internal substructure in multi-parsec 
scale observed clouds (suggested by GMCs' lack of central concentration, and by the weak self- 
gravity of substructures aside from the dense cores and largest clumps [e.g. Bcrtoldi & McKce 
(1992); Williams, de Geus, & Blitz (1994)]), the correspondence between the intermediate-scale 
("clump") structure in real clouds and in our model snapshots may be quite direct. 

Since some ambiguity remains in associating an overall physical length scale with our simu- 
lated models (due to the periodic boundary conditions), we report integrated quantities solely in 
dimensionless units, giving e.g. column densities in units of the mean column density, N = pL. For 
local variables (such as magnetic field strengths) which bear no such ambiguity, we report values 
in dimensionless units and also transform to physical units based on our adopted fiducial density 
and temperature. 






(4) 




In the structural analyses of §§4-6, we report on 
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3. Energy evolution in model clouds 

The early evolution in all the models follows a similar course. Kinetic energy initially decreases 
as the fluid works to deform the magnetic field. The initially non-compressive velocity field is 
transformed into a compressive field, by interactions with the magnetic field and nonlinear coupling 
of the spatial Fourier components. This leads to the development of density-enhanced and density- 
deficient regions, and results in the dissipation of energy in shocks. Fluctuations in the density 
cause fluctuations in the gravitational potential that begin to dominate the dynamics at late times, 
and lead to runaway gravitational collapse for supercritical models. 

To quantify the energetic evolution, we define the kinetic energy 

^K = \j d\ (vl + f 2 + v^,)p, (5) 

the perturbed magnetic energy 

SEB = j^l d^r {Bl + Bl + B^) - Eb,o, (6) 

where Eb,o = L^B^/8 IT is the energy in the mean magnetic field, and the gravitational potential 
energy 

EG = \j d^v{p4>G), (7) 

where is the gravitational potential computed from the Poisson equation modified for periodic 
boundary conditions. 

Figure 1 shows the evolution of Etot = Ek + 5Eb + Eg, Eg, Ek, and Eb (see panels a-d, 
respectively). Prom Fig. lb, it is clear that, consistent with expectations and previous results 
on self-gravitating cloud models with decaying turbulence (Paper III), all but the magnetically 
subcritical (3 = 0.01 model suffers a gravitational runaway. Both of the supercritical models become 
gravitationally bound at time Q.Qtg, corresponding to 6 Myr (nji/j/lO^cm"^) . ^ The 
gravitational runaway time is comparable to that found in lower-dimensional simulations. 

The kinetic energies in all models decay rapidly. After one flow crossing time tf, the kinetic 
energy has been reduced by 73 — 85% compared to the initial value (sec Table 1). The kinetic 
energy is reduced by a factor two after 0.2 — 0.4 flow crossing times (Table 1), with this kinetic 
loss time decreasing toward lower (3 (stronger Bq) because of the more rapid transfer of kinetic 
to perturbed magnetic energy when the Alfven frequency is higher. For GMC parameters (see 
eq. 4), the corresponding dimensional kinetic energy decay time would be 1-2 Myr. The growth 



^The simulations are terminated shortly after the onset of gravitational runaway because the coincident develop- 
ment of low-density regions where va is large causes the Courant-condition-limited timestep to become very short. 
For the /5 = 1 model, the gravitational binding time (the time to reach Etot = 0) is an extrapolation based on the 
evolutions of the /3 = 0.1 and (3 = 1 models shown in Figs. la,b. 
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of magnetic energy stored in these magnetic field fluctuations (due to advection by the turbulent 
velocity field) is apparent in Figure Id; the initial increase is followed by decreasing or flat perturbed 
magnetic energy as the turbulent velocity field decays. The time to reach the maximum perturbed 
magnetic energy lies in the range 0.1 — 0.2 times the Alfven crossing time, similar to what was 
found in lower- dimensional simulations (Paper III). At the point when 5Eb is maximal, it accounts 
for 20-50% of Eturb- The fraction increases with the mean field strength Bq. 

The total "turbulent" energy {Efurb = + ^Eb) secularly decreases in time; after somewhat 
more than half of the initial turbulent energy is lost, the decay approaches a power-law temporal 
behavior with Et^j-b oc (Figure 2). This late-time scaling in non-self-gravitating models of 3D 
MHD turbulence has been noted previously (Mac Low et al 1998; Stone, Ostriker, & Gammie 1998; 
Mac Low 1999). Most of the turbulent losses, however, occur before the onset of this behavior. The 
turbulent decay can be characterized by the time tdec to reduce the turbulent energy by a factor 
two from its initial value. We find that this time is in the range 0.4 — 0.8 flow crossing times (see 
Table 1), comparable to our results from Paper II, and consistent with other findings (Paper III, 
Mac Low et al (1998); Mac Low (1999),Padoan & Nordlund (1999)) that dissipation times vary by 
only a factor ~ 2 over the range of Mach numbers and magnetic field strengths present in GMCs. 
The corresponding dimensional time for turbulent energy decay with CMC parameters is 2-4 Myr. 



4. Density and column density distributions 

A basic statistical property of a real or model cloud is the distribution of density in its con- 
stituent parts. This distribution may be described either by its fractional volume per unit density 
{dV/ dp) or by its fractional mass per unit density {dM/ dp) . Previous analyses of the density distri- 
butions in compressible hydrodynamic turbulence simulations (before gravity becomes important) 
show that when the equation of state is approximately isothermal, the density distribution is close 
to a log-normal (Vazquez- Scmadeni 1994; Padoan, Jones, &: Nordlund 1997; Passot & Vazquez- 
Semadeni 1998; Ostriker, Gammie, & Stone 1999). Scalo et al (1998), Passot & Vazquez- Scmadeni 
(1998), and Nordlund & Padoan (1999) also show that in a medium where the temperature de- 
creases (increases) with increasing density, an extended tail in the density distribution function 
develops at density higher (lower) than the mean density. In the present models we assume an 
isothermal equation of state. This is is a reasonable approximation here since most of the gas is 
contained in condensations at density larger than the mean value where the temperature likely 
varies by less than a factor ~ 2 (see Paper III, also Scalo et al (1998)). 

As described in §1, a key question is whether it is possible to discriminate the magnetic field 
strength in a cloud from its structural properties. Using our present models, we can test how the 
strength of the mean magnetic field affects the observable density and column density statistics. 
For these tests, we choose sets of model "snapshots" from the three decay models in which the 
Mach number M. = (v^/Cg) (or kinetic energy) matches in the three models; because the energy 
evolves at somewhat different rates in the runs with different (3, these times of the snapshots vary. 
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The sets of model snapshots have M. ^ 9,7, 5. 

Figure 3 shows an example of the distributions of volume and mass as a function of volume 
density, p, where p is measured in units of the mean density p = M/L^ in the simulation cube. 
The volume density distributions arc well approximated by log-normal functions, i.e., volume and 
mass distributions in y = log (p/p) of the form 

fv,M (y) = -^L= exp [- (y ± IpI ) 2/2^2] , (8) 

where the upper/lower sign on the subscript applies to the volume/mass distribution; fv,Mdy is the 
fraction of the volume or mass with y in the interval (y, y + dy). It is straightforward to show that 
the mean p and dispersion a are related by \p[y]\ = ln(10)cr[y]^/2 for a log-normal distribution, 
so a[y] = Q.93^/J^[y\.'^ Table 2 (cols. 5,6) gives the values of i^tviv] (= {y)v) and PmIv] (= {y)M) 
for three sets of models at different Mach numbers (where the subscript on the angle brackets 
denotes weighting by volume or mass). In all cases, pv[y] ~ ~A*M[y], consistent with a lognormal 
distribution. 

For a log-normal distribution, the weighted mean and dispersion of the density itself are related 
to the mean of the logarithmic density contrast p[y] using 

log/^\ =2/log^\ ^my]\, (9) 

\p/m \ PI M 



and 




M 



P 

P/ M 



P 

P/ M 



102|M[y]|^ (10) 



106|M2/]I^ (11) 



where "V" and "M" subscripts denote weighting by volume and mass, and (p)y = p. From Table 
2, p[y] is in the range ~ 0.2 — 0.4 for = 5 — 9, implying from equation (9) that the typical mass 
element has been compressed by a factor (p/p)m ~ 2.5 — 6 compared to its unperturbed initial 

value. Because of the log-normal form of the distribution, two-thirds of the matter is within a 
factor lO^-^^Vlt^M (~ 2.7-3.8) above or below the value IQl^'I^'lIp ~ (1.6 - 2.5)p, and 95% is within 
a factor 10^-'^^\/l^[2'll (r^ 7.2 — 14) above or below this value. The volume- weighted rms standard 
deviation in p/p is (lO^^t^l - 1)^/^ 13 _ 2.2), and the mass-weighted rms standard deviation in 
p/p is 102^^[2'1(102'^[2'1 - 1)^2 3 _ 13)_ 

The above results on density contrast may be compared with previous work. In Paper III, 
we found that for 2.5 dimensional models of decaying turbulence with (3 = 0.01,0.1, 1.0, the mean 



^Elsewhere, distributions axe sometimes given as a function of x = ln(p/p); in that case, \tJ,[x]\ = a[x]'^/2 since 
H[y] = //[a;]/ln(10) and a[y] = cr[a;]/ln(10). 
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Table 1. Comparative Energy Evolution 



model 


/3 


EK{tf) /EK,init^ 


5EB/EK{tfT 








B 


0.01 


0.27 


0.51 


0.76 


0.21 


> 0.9 


C 


0.1 


0.15 


1.13 


0.55 


0.31 


0.6 


D 


1.0 


0.17 


0.35 


0.41 


0.37 


0.6 



•^The flow crossing time tf = L/(2£;K,i„it)^/^ ^ 0.07is 

^tdec {tdec) time to reduce the initial energy (kinetic energy) by 50% 

"tbind is the time at which Ex + SEb +EG = 0;tg = (Tr/Gp)^/^ 0.331^ 



Table 2. Comparative Density and Column Density 



snapshot 


/3 


t/ts 


M 












Bl 


0.01 


0.03 


8.8 


-0.25 


0.30 


0.024 


0.033 


0.029 


CI 


0.1 


0.03 


8.8 


-0.27 


0.28 


0.030 


0.037 


0.047 


Dl 


1.0 


0.03 


9.4 


-0.42 


0.37 


0.044 


0.047 


0.061 


B2 


0.01 


0.07 


7.4 


-0.38 


0.38 


0.038 


0.037 


0.046 


C2 


0.1 


0.04 


7.6 


-0.27 


0.29 


0.034 


0.046 


0.048 


D2 


1.0 


0.05 


7.2 


-0.37 


0.34 


0.060 


0.054 


0.065 


B3 


0.01 


0.19 


4.9 


-0.23 


0.21 


0.015 


0.028 


0.021 


C3 


0.1 


0.09 


4.9 


-0.35 


0.37 


0.047 


0.050 


0.056 


D3 


1.0 


0.09 


4.9 


-0.31 


0.33 


0.048 


0.057 


0.063 



'^Volume-weighted or mass-weighted average of the logarithmic density contrast, 

y = \og{p/p): expected sampling error is ^ 10~^. 

^ Mass- weighted average of logarithmic column density contrast, Y = \og{N/pL), 
for projection along x, y, or z; expected sampling error is ~ 10~^ 
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logarithmic density contrast |/i[y]| = 0.2 — 0.5 for Mach numbers in the range 5 — 10, with a weak 
trend toward an increase in the contrast with increasing Mach number, and the largest contrast 
in the strong-field (/? = 0.01) group. For the quasi-steady forced-turbulence models with ~ 5 
reported on in Paper II, the mean logarithmic density contrasts are in the range 0.20-0.28, 

increasing from P = 1.0 to 0.01. Thus, overall, we find comparable values of the density contrast 
in all our analyses of turbulence in those stages where self-gravity is not important. 

Nordlund &; Padoan (1999) and Padoan, Jones, Sz Nordlund (1997) report findings imply- 
ing that, for 3D unmagnetized forced turbulence, fi[y] is related to the Mach number M by 
\l^[y]\ = (l/2)log(l -I- 0.25M'^). For the range of Mach numbers (X ~ 5 - 9) in our Table 2, 
the corresponding values of | would be ~ 0.4 — 0.7, somewhat larger than those we found; how- 
ever, Nordlund & Padoan (1999) remark that they find lower density contrasts when the magnetic 
field is nonzero, which would yield better agreement with our results for magnetized turbulence. 
Analysis of simulations of compressible, isothermal, unmagnetized turbulence in one dimension by 
Passot & Vazquez-Semadeni (1998) suggest a linear rather than logarithmic scaling for fi[y] with 
Mach number A4 and much larger values of the contrast than those found in 3D simulations. This 
may be due to the purely compressive velocity field in ID. 

The present analysis suggests that for non-steady magnetized turbulence, there is no one-to-one 
relationship between the density contrast and the sonic Mach number or other simple characteristic 
of the fiow. There does, however, appear to be a secular increase in the minimum value of the 
contrast with the effective Mach number for magnetized flow, the fast magnetosonic Mach number 
Mf defined by M'j^ = {v^) / {v\ + c^). In Figure 4, we plot the logarithmic contrast factors 
against the value of log(l -I- Alp); the lower envelope of the contrast is found to be fit by = 
0.2(log(l -I- A^p) + 1) for the models studied. There is no similar secular relationship between the 
density contrast and the ordinary sonic Mach number M.. 

Applying similar reasoning to the argument of Passot &; Vazquez-Semadeni (1998), the weak 
relation between the effective Mach number ^AF and the density contrast may be understood 
heuristically as follows. From equations (9) and (11), we may write the mean logarithmic contrast 
in terms of the mass-weighted dispersion in density amplitude as 




In strong, unmagnetized, isothermal shocks, which would occur for fiow parallel to the field, 
the preshock and postshock densities p\ and p2 have {p-z/ Pi) — 1 ~ M.^- For strong isothermal 
shocks magnetized parallel to the shock front and /3 ^ 1, (p2/pi) — 1 is linear rather than quadratic 
in Aip^ approaching ^pivjvA ~ \f2J^p. If the typical shock jump compression factor determines 
the rms dispersion in the density, then the term in square brackets in equation (12) would scale 
between quadratically and quartically in A4f for a range of /3 and shock geometries (noting that 
M.F ^ M. ioY P large). The real situation is of course more complicated. It is interesting, however, 
that the slope ~ 0.2 of the lower envelope of the < log(/9/p) > vs. log(l -|- A4p) relation does fall 
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in the range between 0.17 and 0.33 suggested by this heuristic argument. The fact that this lower 

envelope lies closer to the (shallower) slope corresponding to parallel-magnetized shocks indicates 
that the model turbulent clouds do not invariably evolve to be dominated by (more compressive) 
flows aligned with the mean magnetic field. 

Because of the potential for direct comparison with observation, it useful to examine the 
distributions of column density N. In particular, we would like to ascertain if the distribution 
depends on the mean magnetization. The distribution of column densities can be described by 
the fractional area, dA/dN{s), or fractional mass, dM/dN{s) per unit column density, where s is 
the orientation of the line of sight through the cloud. In Figure 5, we compare the distributions 
of projected area and mass as a function of column density for model snapshots (B2,C2,D2 from 
Table 2) with matched Mach numbers and different values of the mean magnetic field strength. 
Although the statistics are poorer than for the distributions of volume density, the column density 
distributions are also approximately log-normal in shape. Thus, the column density distributions 
can be described using the same form as equation (8), but replacing y ^ Y = log{N/N), where 
N = pL. The mean and dispersion may depend on the projection direction s, so ^[y\ and 
^\y\ ~^ <^s\Y]. To the extent that the distributions follow log-normal forms, the formal relations 
(9)-(ll) would apply, with p/p — > N/N and area-weighting replacing volume-weighting. 

In Table 2, we list the values of the mass-weighted mean of the logarithmic column density 
contrast (|;Um;s[^]| = (log(A^/^))j^,^.^) for the different model snapshots in each of three projection 
orientations s = x,y, z. From the data in the Table, the projections in the various directions for 
any model snapshot yield somewhat different statistics (mostly 10-20 % differences in a[Y] and 
twice that in ju[y]); the projection along the magnetic field tends to give slightly lower contrast 
than the two perpendicular projections. Differences between the two perpendicular directions {z 
and y) arc simply a result of specific realizations of random initial conditions. Because the models 
have the same initial velocity perturbation realization, they will have similar evolved structure to 
the extent that the magnetic fields only weakly affect the dynamics - this explains, for example, 
why models C and D both have larger contrast for z projections than y projections. 

Notice that models with the strongest magnetic field tend to have lower column density con- 
trasts than models with the same Mach number and weaker mean Bq (20-50% differences in a[Y] 
for most sets). This effect is most pronounced for the Mach-5 set (B3, C3, D3; this set has a 
factor two [three] difference in the a[Y] [iJ,[Y]]). Padoan &; Nordlund (1999) previously pointed out 
that column density contrasts may be larger in weaker mean-Bo models. Our results confirm this 
tendency, although we find that the effect is relatively weak in magnitude, and does not hold in all 
cases (see e.g. the results for snapshots B2 and C2 in the Table). 

Overall, the range of mean logarithmic column density contrasts in Table 2 is fJ-[Y] ~ 0.015 — 

0. 065, corresponding to typical mass-averaged column density in the range {N/N)m = 1.07 — 1.35, 

1. e., only a modest enhancement over the average in a uniform cloud. The range of logarithmic 
column density contrasts is an order of magnitude lower than the range of mean logarithmic density 
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contrasts. This is understandable, since each column contributing to the distribution samples a 
large number of over- and under- densities along the line of sight. The column density distributions 
still require a density correlation length over a significant fraction of the box size L along the line 
of sight, however; otherwise the dispersion in column densities would be wiped out by line-of-sight 
averaging. 

This can be seen more quantitatively as follows. Each of ha columns that contributes to the 
distribution is created by taking the sum of densities in Ug cells along the line of sight. Prom 
the Central Limit Theorem, we know that if the density in each cell along the line-of-sight were 
an independent random variable, then for Ug large, the distribution of column densities would 
approach a Gaussian - rather than log-normal - shape, with (area- weighted) mean of N/ N equal 

- _ _ _ _ _ —1/2 

to 1 (where N = pL), and (area- weighted) standard deviation in N/N equal to Us times the 
(volume-weighted) standard deviation in p/p. For a log-normal volume density distribution obeying 
equation (10), an assumption of independent sampling along the line-of-sight would therefore predict 
an (area- weighted) rms deviation of N/ N from unity given by 

^ . _ _ . 1/2 



Gauss _ ^iQ2\tj.[y]\ _ ly ^ ^i^^ 



'N/N 



with typical sampling error ~ ^ NjN I determining the mean and dispersion of N/N. For 
\p[y] \ ~ 0.2 — 0.4, the expected standard deviation in N/N would be ~ 0.08 — 0.14, with sampling 
error ~ 0.0003 — 0.0006, if the line-of-sight cells were all independent. In fact, using the area- 
weighted equivalent of relation (10) for the log-normal {not Gaussian) column density distribution 
that is evidently produced, the area- weighted standard deviation in N/N is 

^lo9_norm ^ f^^Q2HY]\ _ ^2 ^ ^^^^ 

or approximately y^21n(10)/x[y] for fM[Y] « 1. For our tabulated values, this is in the range 
0.27 — 0.59, significantly larger (by hundreds of times the sampling error) than would be predicted 
by assuming uncorrelated values of the density along any given line of sight. Thus, both the non- 
Gaussian shape and the breadth of the dispersion of the column density distributions argues that 
the volume densities are not independent but are correlated along any line of sight - as indeed 
should be expected since there are large coherent regions of density created by the dynamical flow. 

We speculate that it may be possible to understand the dynamical process behind the de- 
velopment of the log-normal column density distribution following similar reasoning to the argu- 
ment of Passot & Vazquez-Semadeni (1998) for the development of a log-normal volume density 
distribiition. They argue that if consecutive local density enhancements and decrements occur 
with independent multiplicative factors due to independent consecutive velocity compressions and 
rarefactions, then the log of the density in some position is the sum of logs of independent en- 
hancement/decrement factors; this would yield a lognormal density distribution if there are many 
independent compressions/rarefactions, each sampling independently from the same distribution of 
enhancment / decrement factors . 
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Suppose, similarly, that the gas along any line of sight is subject to multiple independent 
compression/rarefaction events; since the compression/rarefaction axes are not in general along the 
line-of-sight, column density on a given linc-of-sight is not conserved. Each compression/rarefaction 
event which produces a local change in the volume density by a factor X affects only a fraction / 
of the column of gas, resulting in an effective enhancement /decrement factor for the column closer 
to unity than X. A simple model would be to suppose that each event i independently produces 
a change in the column density by a factor X'^ = (1 — fi) + f^Xi = Xi — {1 — fi){Xi — 1) (taking 
the fraction 1 — /j of the gas in the column at unchanged volume density and the fraction fi at 
volume density enhanced/decreased by a factor Xj). If > 1 (respectively, Xi < 1), then X'- < Xi 
(respectively X- > Xi). The logarithm of the column density contrast would then be a sum of terms 
logX'-; taking these as random variables, the resulting distribution would be log-normal (assuming 
a large number of [spatially overlapping] successive events). Since each X'^ is closer to unity than 
Xi, the mean and dispersion of the logarithmic column density distribution arc expected to be 
smaller than those of the volume density distribution. Although it would be interesting to test in 
detail whether this sort of heuristic model could be refined and used to relate projected density 
distributions to volume density distributions, the potential for finding a unique inversion (even in 
a statistical sense) is limited by the many degrees of freedom associated with the extended spatial 
power spectrum producing the compressions. 

For the power-law input turbulent spectrum that we adopt, the spatial correlations which 
produce the column density distribution occur at sufficiently large scale that the distributions are 
not, except at columns N much larger than the mean, very sensitive to the "observer's" resolution. 
For example. Figure 6 shows the statistics of column density for one model at the full resolution of 
the simulation, and for resolving power reduced by factor of four by averaging the column density 
values within squares of edge size four times that of simulation cells (so that each "pixel" has 
sixteen times the area of a projected simulation cell). The overall shapes and mean values of the 
distributions are quite comparable. At column densities much larger than the mean, of course, 
the distributions become sensitive to resolution because of the scarcity of regions with the highest 
column density; averaging these with their lower-column-density neighbors results in a cutoff of 
the distribution at lower N. A related point for observed ^^CO data was discussed by Blitz & 
Williams (1997). They showed that the distribution of the number of cells in position- velocity 
space as a function of TA/TA,max in the cell flattens as the linear resolution scale increases, due to 
the smearing-out of the highest-column regions. We have verified that the distribution of number 
of projected cells with N/Nmax similarly becomes flatter if the map of projected density is averaged 
over grids with increasing cell size. 

Because the periodic boundary conditions introduce an effective correlation in the density along 
the line-of-sight at scale L, a potential concern might be that the typical column density contrast 
might be enhanced by introducing "artificial" coherence along lines of sight. ^ To investigate this 



^We thank E. Vazquez-Semadeni for noting this point. 
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effect, we have evaluated sets of "half-column" density distributions by summing only over distances 
L/2 along the line-of-sight. In general, the resultant half-column distributions are still lognormal 
in form (although noisier), with larger means and dispersions than those found for the full-column 
integrations. Figure 7 shows one such set of distributions, obtained from the z > L/2 "front half" 
of the volume snapshots B2, C2, D2. This result suggests that the coherent volume density regions 
responsible for the lognormal column density distribution in fact have intermediate scale - they 
are much larger than the cell size, but significantly smaller than the overall size of the box. In 
this situation, one would expect that a factor two decrease in the number of (multi-cell) correlated 
regions along the line-of-sight would produce a factor \/2 increase in ctj^^j^, corresponding to a 
factor ~ 2 increase in (cf. eq. 14). Indeed, we find the half-column values of are 

typically larger than the full-column values by a factor ~ 1.5 — 2, supporting this interpretation. 

The robustness of the column density distribution to resolution changes makes it a viable 
statistic for comparing simulations to the observable properties of turbulent clouds. Such compar- 
isons are a test of the idea that much of the moderate-density "clumpy" structure in molecular 
clouds may be produced by turbulent stresses. Preliminary results are promising; for example, we 
have compared the distribution of the extinction data values from the dark cloud IC5146 (Lada, 
Alves, &: Lada 1999) with column density distributions from our simulation snapshots. Figure 8 
shows that the cumulative distributions are indeed remarkably similar in form (although this par- 
ticular real cloud has a slightly larger dispersion than our models have). Unfortunately, however, 
the column density distribution is determined by more than just a few simple global parameters. In 
some circumstances, there may be as much variation in the column density distributions between 
the same cloud viewed at different orientations as there is in two clouds with the same turbulent 
Mach number but a factor ten difference in the mean magnetic field. This large "cosmic variance", 
and the relatively weak variation with parameters of and a compared to their scatter, make it 
unlikely that it will be possible to estimate individual clouds' mean magnetic field strengths, for 
example, from column density distributions alone. 



5. Linewidth-size relations and projection effects 

An important way of characterizing the kinetic structure in turbulent clouds is to measure the 
distribution of the velocity dispersion vs. the physical size or mass of the regions over which it 
is averaged. Means over these distributions then represent "linewidth-size" relations. The regions 
over which velocity dispersions are averaged in observed clouds are often apparent "clumps". At 
the most basic level, an apparent "clump" in a cloud or projected cloud is a spatially connected, 
compact region that stands out against the surrounding background. In any hierarchical structure, 
clumps will contain other smaller clumps, and in general the identification of clumps is a resolution- 
dependent procedure. Starting from the fundamental concept of a "clump" as a region of contrast 
(ROC) on a given spatial scale, we have developed a simple algorithm to identify and characterize 
the ensemble of projected ROCs at multiple scales, so as to explore the scaling of kinetic properties 
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with physical size.^ 

The procedure is as follows: First, we choose a size scale s (here, a factor 2" times the simulation 
grid scale, where n = 2 — 8). We then divide the projected cloud into zones of area s^. Within each 
zone, we compute the mean projected surface density as the total zone mass divided by s^; we also 
compute the mass- weighted mean surface density S[s] for the set of zones on scale s. We label a 
zone as a ROC on scale s if its surface density is at least a factor fc times 5][s]. Typically we use 
/c = 1, but the results are not qualitatively sensitive to this choice; we note that (i) regions above 
the mean column density at a given scale occupy less than half the area due to mass conservation, 
and (ii) since S[s] increases with decreasing s, the ROCs on a given scale would appear "by eye" 
to stand out against the background even with fc = 1. For each projected ROC, we also compute 
the (mass- weighted) dispersion of the line-of-sight velocity (t„; this represents the "linewidth" for 
a region of projected area s^. 

We are now in a position to examine the correlations among linewidth a^, mass M, and spatial 
size s for our ROC collections. In a data set based on molecular line emission, the contributions 
from any local region would depend on the local excitation rather than simply being proportional to 
the amount of matter present. For the analysis described below (except as noted), we only include 
contributions from material if its local density (mass/volume) is at least equal to pmin = 3/0, as 
a simple way of selecting material in the range of densities that contribute to common molecular 
lines. ^ 

Figure (9) shows an example of how the ROCs at multiple scales are distributed on the map 
of model snapshot B2 projected in the z direction (Figure 22 shows a colorscale image of the 
column density for the same snapshot projection). For the ROC ensemble shown in the figure, we 
compute masses, velocity dispersions, and values of the so-called virial parameter a = 5a^s/GM 
Bertoldi & McKee (1992). In Figure (10), we plot the values as a function of (linear) size scale 
and/or mass. We also evaluate least-squares linear fits to dlog((7„)/dlog(s), d\og{M) / dlog{s) , 
d\og{ay)/dlog{M), and dlog(Q;)/(ilog(s); the respective values in this example are 0.09, 1.87, 0.06, 
and -0.4. 

From parts (a) and (c) of Figure (10), it is clear that although the there is a mean increase 
in velocity dispersion with mass and linear size, there is a great deal of scatter as well. The upper 
envelopes of the velocity dispersion distributions in fact even decrease as a function of increasing 
M and s; the lower envelopes increase more steeply. The distribution of a vs. M also shows large 
dispersion, with a nearly-flat lower envelope and an upper envelope showing a decrease in a with 
M. The M vs. s distribution has a relatively low dispersion. 



^In Gammie, Stone, & Ostriker (2000), we use an alternative approach to define clumps and characterize their 

properties. 

'^Realistically, of course, the contribution to observed lines depends on more than the local density; due to radiative 
transfer effects, it might even be possible for lower-density material to contribute more efficiently than higher-density 
material if its emission occurs in line wings and suffers less absorption. 



-19- 



Many of the features evident in Figure (10) can be understood by reference to the scahng 
properties of the underlying three-dimensional distribution, together with the effects of projection 
onto a plane. First consider the log((J^,) - log(s) distribution. The procedure we have used to identify 
ROCs in the projected plane also can be used in the 3D data cube itself; we can then compute 
the mass and velocity dispersion for each 3D cell of edge size s that meets the contrast criterion. 
In Figure (10 a), we show how the mean velocity dispersion for these 3D cells depends on size 
scale. Interestingly, this curve traces fairly closely the lower envelope of the distribution of o"„ vs. 
projected size for ROCs on the projected plane. Thus, for nearly all projected regions of area s x s, 
the majority the velocity dispersion can be attributed to the superposition along the line-of-sight 
of many regions of volume s x s x s with different mean velocities. The relatively weak dependence 
of mean linewidth on projected size (or mass) simply reflects the ubiquity of "contamination" by 
foreground and background material. Previously, Issa, MacLaren, & Wolfendale (1990) and Adler 
Sz Roberts (1992) have made a related point that inferred broad linewidths of apparently quite 
massive GMCs may arise due to overlapping in velocity space of narrower velocity distributions 
from individual smaller clouds superimposed along the line of sight. Relatively steep increases 
of linewidth with size, as reported by Larson (1981) and subsequent authors, may be obtained 
in observations provided that a structure is distinguishable from its surrounding by a sufficient 
density or chemical contrast; these steeper laws correspond to what we measure with our 3D ROC 
procedure (solid lines in 10a, 11a). 

In Fig. (10b), we show the distribution of mass with projected size for the ROCs; the mean 
logarithmic slope is nearly equal to two - rather than three, as would be the case for compact 
objects with three comparable dimensions. The mean slope is close to two simply because each 
ROC samples along the entire line-of-sight so that mass is nearly proportional to projected area; 
note, however, that at small scales, the masses can lie considerably above the mean fit. 

It is interesting to compare the virial parameter a vs. M distribution shown in Figure (10) 
with the analogous plot presented by Bertoldi & McKee (1992) analyzing the properties of apparent 
clumps in four different observed clouds. For the data sets considered in that work, linear fits to 
the logo; vs. logM relation gave slopes between —0.5 and —0.68. For the model data shown in 
Fig. (10) (and for our other snapshots as well), the mean fit has a somewhat shallower slope. But 
the upper envelope of this distribution (and those for other snapshots) has slope ~ —0.5 to —0.6. 
We can understand this upper envelope as follows: First, the largest velocity dispersions at a given 
projected scale (cf. Fig. 10a) are nearly independent of scale (typical logarithmic slope is ~ to 
—0.1). With this, together with the mass scaling nearly as s^, the result is an upper envelope of 
a oc M~°-^ — M~°-^. The relatively fiat lower envelope of the a vs. M distribution can be explained 
by the projected ROCs that sample the (t.^, oc s"'^ lower-envelope of the linewidth-size distribution 
(following the true 3D linewidth-size relation) , together with the approximate M oc scaling. 

All of the other model snapshots show qualitatively similar distributions of the kinetic parame- 
ters for ROCs to those shown in Figure (10). For example, we show the same distributions obtained 
for a weak-magnetic-field model snapshot (D2) in Figure (11); qualitatively, all of the kinetic seal- 
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ings are quite comparable to those obtained for the strong-magnetic-field model. In general, for the 
model snapshots in Table 2, the projections parallel to the magnetic field axes yield slightly stronger 
increase of linewidth with size than do the other projections. For projections perpendicular to the 
mean field, the ranges in the fits for the different snapshots are dlog(c7„)/dlog(s) = 0.07 — 0.12, 
dlog(cr^)/cilog(M) = 0.03 - 0.08, and dlog(a)/(ilog(s) = -0.45 to -0.34 (using the same mini- 
mum surface density contrast factor /c = 1 and Pmin = For projections parallel to the mean 
field, the respective ranges for these fits are 0.11 - 0.19, 0.06 - 0.12, and -0.40 to -0.29. The 
fits to (ilog(M)/dlog(s) have a very small range, 1.83 — 1.89, for all projections (using /c = 1 and 

Pmin — 

The results depend weakly on the definition of a ROC, and in particular on Pmin- Reducing 
Pmin tends to yield flatter slopes for d\og{a.u) / d\og{s) and dlog(o-^)/(ilog(M) (because velocity is 
anticorrelated with density, so that additional low-density material along the line of sight increases 
the dispersion closer to the maximum), and steeper slopes for d\og{M) / d\og{s) (approaching 2, 
the limiting form for uniform column density), and for c/log(Q;)/dlog(s) (approaching —0.5, the 
limiting form for velocity dispersion independent of size and uniform column density). Increasing 
Pmin has the Opposite effect. The changes in slopes come about mainly from variations in the loci 
of the lower envelopes of the distributions when Pmin varies; the upper envelopes change very little, 
since they reflect the kinetic properties of ROCs which sample through the largest possible portion 
of the model cloud. 

Because the projected ROC identification algorithm does not take into account any line-of- 
sight information for the material in any projected region, it should not be surprising that the 
velocity dispersions for projected regions can be much larger than the velocity dispersions for 3D 
cubes with the same projected size. One might argue that foreground and background material 
extraneous to a principal condensation could easily be removed based on velocity information, so 
that structures identified as contrasting regions in observed molecular line I — h — v data cubes 
would truly represent spatially coherent structures. Examination of the line-of-sight velocity and 
line-of-sight position distributions for individual projected ROCs, however, suggests that it may in 
fact be difficult to eliminate foreground/background contamination. 

To illustrate the problem. Figure (12) shows histograms of line-of-sight velocity (equivalent to a 
line profile for an optically-thin tracer uniformly excited above Pmin = 3p) for regions of projected 
linear scale s = L/8. Although some of the line shapes are irregular, none of those meeting 
the ROC criterion (in this example) are clearly multicomponent distributions. For comparison, in 
Figure (13), we show the distribution of mass with position along the line of sight. Evidently almost 
every region - both ROCs and non-ROCs - has multiple spatial components along the line-of-sight. 
Figures (14) and (15) show the same distributions for spatial regions at higher resolution; again, 
almost all velocity profiles are single component, while spatial distributions are multicomponent. 
By dividing our data cubes in half and computing velocity histograms separately for the "front" 
and "back" halves, we have checked that the ubiquity of single-component velocity distributions 
is not an artifact of periodic boundary conditions. We have also checked that the phenomenon of 
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multi-spatial component/ single- velocity component ROCs is still prevalent even when the density 

threshold pmin is set higher; for example, Figure 16 a,b shows the velocity and position distributions 
for the same model as before, but now with pmin = 10/?. A complementary phenomenon that we 
have also identified in several model projections is that multiple velocity components in a given 
ROC may correspond to a single extended spatial component - as a consequence, for example of 
two "colliding" clumps being viewed during a merger along the line-of-sight. 

The general lack of correspondence between structures in position and velocity space in ISM 
models has previously been noted, based on various sorts of analyses. For example, Adler & 
Roberts (1992) analyze the model galactic disks generated from two-dimensional N-body cloud- 
fluid simulations, and show that apparent single "clouds" in longitude-velocity space are often 
highly extended along the line-of-sight, and that what appears to be a single GMC in a spatial 
plot may be assigned to multiple "clouds" in longitude-velocity space. Pichardo et al (2000) show 
that the morphology of structures in position-position-velocity space (equivalent to channel maps) 
in their 3D MHD simulations is more strongly correlated with velocity structures in physical space 
than with density structures in physical space. 

The current analyses and previous work on this question do not treat molecular excitation 
and radiative transfer in detail. Studies that do include these complex effects will be required 
to reach definitive conclusions on the relation between maps of molecular lines and 3D physical 
density-temperature-velocity cubes. If spatially compact regions have substantially higher molec- 
ular excitation than more diffuse surroundings due to line trapping, then it is still possible that 
velocity information could be used to separate spatially-connected clumps from foreground and 
background material. Large amplitude rotation of clouds, if present, would also help to differentiate 
superposed line-of-sight clumps in the velocity domain. Potentially, methods that use specific infor- 
mation about spectral line shapes (e.g. Roslowsky et al (1999)) may also be adapted to discriminate 
spatially-separated regions. The present simplified analysis suggests, however, that foreground and 
background material may at least significantly increase the dispersion in the linewidth-size distri- 
butions for clumps identified from molecular emission data cubes (e.g. Williams, de Geus, k. Blitz 
(1994), Stutzki & Glisten (1990)). 



6. Magnetic field distributions and simulated polarization 

We now turn to magnetic field structure, and begin by considering how the distribution of 
the magnetic field varies for models with different mean magnetization. As seen in Figure Id, 
the rms magnetic field strength initially increases, due to the generation of perturbed field by 
velocity shear and compression. The distribution of the individual components of B for matched 
Mach number model snapshots B2,C2,D2 is shown in Figure 17. The dimensionless field strength 
B = 5/(47rpc^)^/^ that we report can be converted to a physical value using 



B = 1.4 X i3 X (T/10K)V2(„^Jio2cm-3)V2 ^q. 



(15) 
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As illustrated by the figure, the component distributions are more nearly Gaussian for the case 
of stronger magnetic fields; this is true for all of the model snapshots, although the distributions 
in the high-/? (low-i?o) cases do become more Gaussian in time. For the weak-field models, the 
dispersion in each component of the magnetic field is larger than the mean field component. 

Because magnetic fields are measured via the Zeeman effect with different atomic and molecular 
tracers in different density regimes, it is interesting to analyze how the mean field strength in 
simulations may depend on density. Since the magnetic field is weaker and less able to resist being 
pushed around by the matter in the /3 = 0.1, 1 (C and D) simulations, one expects that the field 
strength will have stronger density dependence for these models than for the /3 = 0.01 simulation. 
This is indeed the case, as can be seen in Figure 18. Particularly at densities below the mean, the 
magnetic field strengths in the high-/? models are strongly density dependent; the low-density slope 
of dlogB/dlogp for these models is near the value 2/3 associated with a constant ratio of mass to 
magnetic flux and isotropic volume changes. 

At high densities (above ~ lOp) the relatively flat slope of the P = 0.01 model increases, 
becoming comparable to the slopes of the /3 = 0.1, 1 models. Figure 19 shows the high-density B vs. 
p dependence, for various model snapshots; fits for fiducial density in the range 10^ — lO^cm"^ 
(i.e. p/p= 10 to 100) yield slopes 0.3-0.7 for dlogB/dlogp. Only the Mach-9 P = l model yields 
a high-density-regime slope as steep as the isotropic contraction limit. The other snapshots have 
slopes 0.3-0.5, which may be compared with the slope 0.47 found from a compilation of Zeeman 
measurements at high densities nj^j = 10^ — 10 ''cm"'^ (Grutcher 1999). The values of the mean 
in any density regime generally increase with increasing mean net magnetic flux Bq (i.e. decreasing 
/?), but because there is significant dispersion about the mean B^, there is considerable overlap of 
the 1 — cr deviation regions among the different model snapshots (Fig. 18) . 

The numerical results on the B vs. p relation presented by Padoan &; Nordlund (1999) (see 
their Fig. 7) are qualitatively similar to our results, with some differences apparent at the high 
density end. The lower Mach number in their low-/? model compared to their high-/? model likely 
accounts for its relatively weaker increase of B with p at high density, compared to our results. We 
also differ with those authors regarding the astronomical implications of the numerical results. In 
particular, we do not attempt a comparison of the low-density end of the B vs. p distributions with 
observations made in the diffuse ISM, because (a) the physical regime modeled by the simulations 
is not appropriate for the diffuse ISM (where thermal pressure is comparable to, rather than much 
smaller than, {pv'^)); and (b) the transformation from simulation to physical variables for local 
magnetic field values involves multiplying by the mean magnetic field Bq on the largest scale, and 
this need not be the same in the diffuse and cold ISM (/? parameterizes this mean field strength) . We 
conclude that the B vs. p relations obtained from simulations do not at present constrain the value 
of /?. At high densities, all models (either weak or strong Bq on the large scale) yield slopes which 
are consistent with high-density molecular Zeeman observations. At very low densities, where the 
predictions of models with varying (3 do differ, estimating B within clouds would be difficult, since 
HI Zeeman observations probe the high columns of foreground and background material, rather 
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than the low column of cloud material (although velocity information may help; cf. Goodman &: 
Heiles (1994)). The field strengths in the low density regions within molecular clouds may in fact 
be systematically higher than those at comparable density in the diffuse ISM. 

For all the snapshots, there is significant dispersion in the total magnetic field strength. In 
addition to this overall dispersion in magnitude, there is a dispersion in the magnetic field vector 
direction which increases with decreasing strength of the mean field component Bq, simply because 
fixed amplitude fluctuations have larger relative amplitudes compared to a weak mean magnetic 
field. The dispersion in field directions has important consequences for any observational measure- 
ment of the mean magnetic field via Zeeman splitting. Observations of Zeeman splitting at any 
position on a map yield the line-of-sight average value for the line-of-sight magnetic field, weighted 
at each point along the line-of-sight by the local excitation. When a given line of sight has many 
fluctuations in the direction of the magnetic field, the average value of (-B^os) will be small, even if 
individual local components of the field have large magnitudes. 

To demonstrate how the averaged line-of-sight field components vary with mean field strength 
and observer orientation, we depict in Figure 20 an overlay of (-B/os) on the column density for 
three model snapshots with matched Mach number. In Figure 21, we plot the values of (Bios) 
vs. column density of dense gas. The Figures show, unsurprisingly, that the line-of-sight-averaged 
magnetic field strengths are greatest when the mean field Bq is largest and is oriented along the 
observer's line-of-sight (top left panel). For the weaker-field models, the average line-of-sight field 
is lower, and there is larger dispersion. For all the snapshots, there is considerable dispersion in the 
values of {Bios) on the map, and the largest values do not correspond to the positions of highest 
column density; in fact, there is some tendency of line-of-sight-averaged field to anticorrelate with 
column density. Thus, although the local field strength \B\ increases with density (cf. Figs 18, 19) 
and may be much larger than the volume-averaged mean field Bq for the entire box, line-of-sight 
superpositions of non-aligned vector components produce average line-of-sight field strengths closer 
to the large-scale volume- averaged value. 

It is well known that it is difficult to detect the Zeeman effect in molecular clouds (e.g. Heiles 
et al. (1993)) because the frequency splitting is small when the field is weak. This, coupled with 
the possibility (cf. Fig. 21) that an impractically large number of measurements might be required 
to obtain statististically-significant results for the large-scale field, underscores the importance 
of supplementing programs of direct detection with other methods for estimating the mean field 
strength. Long before direct Zeeman detections were first made, Chandrasekhar & Fermi (1953) 
estimated mean spiral-arm field strengths Bq from the mean gas density, line-of-sight velocity 
dispersion, and the dispersion in orientations of the magnetic field in the plane of the sky. The 
field line orientation is taken to be traced by the polarization direction for background stars, which 
occurs provided that the dust grains producing the intervening extinction are aligned with short 
axes preferentially parallel to B, and so preferentially extinguish linear polarizations perpendicular 
to B. 
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The Chandrasekhar & Fermi (1953) (hereafter CF) estimate is based on the fact that for hnear- 
amplitudc transverse MHD (Alfven) waves, Bp = ^4TTp\5v\/{\6'B\/Bp). Here Bp is the projection of 
the mean magnetic field on the plane of the sky, and (5B and 6v are the components of the magnetic 
and velocity perturbations in the plane of the sky transverse to Bp. If the interstellar polarization 
is parallel to the local direction of Bp, then the ratio |(5B|/5p in the denominator may be replaced 
by the dispersion Scf) in polarization angles (for small angle/low amplitude perturbations). With 
the further assumption that the true velocity perturbations are isotropic, then the dispersion in 
the transverse velocity |5v| is equal to the rms line-of-sight velocity dvios- We thus obtain 

Bp = ^vios S4>-'^ (16) 

where, for CF's field model, Q = 1. Modifications to the CF formula allowing for inhomogeneity and 
line-of-sight averaging are discussed by Zweibel (1990) and Myers & Goodman (1991), respectively. 
Both of these effects (and others; see Zweibel (1996)) tend to reduce Q. 

Potentially, the CF "polarization-dispersion" method can be used to estimate plane-of-sky 
magnetic field strengths on scales within turbulent interstellar clouds. It may also be possible to 
combine these results with Zeeman measurements to estimate the total magnetic field strength 
(Myers Sz Goodman 1991; Goodman Sz Heiles 1994). In order to evaluate the ability of the CF 
method to measure mean plane-of sky field strengths, we provide a first (simplified) test of it using 
our model turbulent clouds. For this test, we have created simulated polarization maps for each 
cloud by integrating the Stokes parameters along the line-of-sight over a projected grid of positions, 
assuming the polarizability in each volume element is proportional to the local density. The details 
of this procedure, together with a more extensive discussion of simulated polarization distributions, 
will appear in a separate publication (Ostriker et al 2000). 

For two projected model snapshots (B2 with /? = 0.01 and D2 with (3=1 projected along z), 
Figures 22 and 23 show examples of the polarization maps overlaid on color scale column density 
maps. The analogous map (not shown) for the model C2 (/? = 0.1) looks quite similar to Figure 
23. From the figures, it it immediately clear that the model with a stronger mean magnetic field 
Bq has more ordered polarization directions and larger typical values of the fractional polarization, 
compared to the model with a weaker mean magnetic field. These trends are as expected: a 
weaker mean field has lower tensile strength, so that for a given level of kinetic energy the Reynolds 
stresses will produce larger fractional perturbations in the magnetic field - corresponding to larger 
fluctuations in projected line-of-sight averaged position angle. Also, because of the larger dispersion 
in local polarization directions along any line of sight, cases with weaker mean magnetic fields will 
show lower net polarization through the cloud (from the line-of-sight averaging of the varying 
local vector directions). While local (line-of-sight averaged) polarization directions may have any 
orientation with respect to local projected surface density, there is some tendency for the large- 
scale projected density and large-scale polarization directions to align in the high-/? (but not low 
f3) models, because the magnetic field and density are both strongly sheared and compressed by 
the large-scale, large-amplitude velocity field. 
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In Figure 24, we show the distributions of polarization angle at various "observer" orientations 
for models with matched kinetic energy and mean magnetic field strengths at three different levels 
{P = 0.01,0.1,1, corresponding to fiducial B = 14,4.4, and lAiiG from cq. 15). As is clear 
from the Figure, only the strong-field model has significantly correlated directions in the simulated 
polarization vectors. This is expected, since only this model has perturbed magnetic energy lower 
than the mean magnetic energy; the ratios are SB^/Bq are 0.27, 4.0, and 12, respectively, for the 
snapshots presented. 

For the cases shown in Figure 24 where the angle dispersion d(f) is 25° or less (i.e. the (3 = 0.01 
projections at i = 0,30,45, and 60°), we have compared the known value of the mean plane- 
of-sky magnetic field Bp = Bq cos(i) with the Chandrasekhar-Fermi estimate. We find that Q 
(see equation (16)) is in the range 0.46-0.51. This suggests that the CF estimate, modified by a 
multiplicative factor ~ 0.5 to account for a more complex magnetic field and density structure, can 
indeed provide an accurate measurement of the plane-of-sky magnetic field when the polarization 
angle fluctuations are relatively small. The method fails, however, when the polarization angle 
fluctuations are large. We will present a more comprehensive analysis of this promising diagnostic 
in a separate publication. 

7. Summary and discussion of structural analyses 

With modern high-performance computational tools, it is possible to create and evolve simu- 
lated dynamical representations of turbulent, magnetized clouds at comparable plane-of-sky spatial 
resolution to that of radio-wavelength observational maps of GMCs. This paper reports on the 
properties of a set of such simulations. 

We start by briefly summarizing (§3) the results on energy evolution in our simulations. We 
conflrm the conclusions from our previous work that turbulent decay is rapid even in magnetized 
models, finding that an interval of only 0.4-0.8 flow crossing times is sufficient to reduce the total 
turbulent energy by a factor two from its initial value; the corresponding physical time for GMC 

parameters is only a 2-4 million years. We also confirm that in situations where turbulence is not 
replenished, the criterion for a cloud to collapse gravitationally depends only on whether it is sub- 
or super-critical with respect to its mean magnetic field; the characteristic collapse time in the 
latter case is ~ 6 Myr for GMC parameters. 

Following the presentation of energetics, the bulk of the paper (§§4-6) is concerned with de- 
veloping tools for structural analyses and applying them to our simulated data cubes. Although 
simplified in their treatment of small scales (ambipolar diffusion is neglected) and thermal proper- 
ties (a constant gas temperature is assumed), the 3D data cube "snapshots" from our numerical 
experiments provide a detailed portrait of the density, velocity, and magnetic field structure in the 
simulated clouds. This structural portrait is dynamically self-consistent in that it is an instanta- 
neous solution to the full time-dependent MHD equations: the density and magnetic field variables 
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have evolved in response to a (time-dependent) turbulent velocity field, which itself has evolved 
subject to gas pressure gradient forces, magnetic stresses, and self-gravity. 

Model cloud snapshots from simulations provide a unique opportunity to (i) explore the intrin- 
sic character of 3D structure in magnetized gaseous systems subject to supersonic turbulence, and 
(ii) determine which aspects of the observed properties of GMCs (from 2D planc-of-sky integrated 
maps oi I — b — V data cubes) can be explained as a manifestation of their internal turbulence. 
The possibilities for such exploration are enormous; for practical purposes, we have limited the 
scope of this paper to three groups of analyses. We consider: (1) the distributions of mass, volume, 
and area as functions of volume density and column density (§4); (2) the distributions of velocity 
dispersion, mass, and virial parameter a as a function of the spatial scale for zones in projected 
maps and cells in 3D cubes (§5); (3) the distributions of magnetic field strength vs. local volume 
density, line-of-sight-averaged line-of-sight magnetic field vs. column density, and distribution of 
simulated polarization angles (§6). For each of these analyses, we compare sets of cloud snapshots 
in which the turbulent Mach number is matched, and the large scale mean magnetic field strength 
Bq varies by a factor ten, also allowing for different "observer" viewing angles. The rms tubulent 
velocities for the model snapshots are = 1 — 2 kms~-^, and the mean magnetic field strengths are 
Bq = 1.4 — 14 iJ,G, assuming fiducial GMC parameters for volume-averaged density = 100 cm~^ 
and temperature T = lOK. 

The main results of these structural analyses are as follows: 

1. The distribution of volume densities follows an approximately log- normal form, with densities 
of typical mass elements {p)m compressed by a factor ~ 3 — 6 times the volume- aver aged density 
p = M/L^ for our sets of snapshots with Mach number M. = ay/cs in the range 5 — 9 (see Fig. 3 
and Table 2). This typical density contrast is comparable to that inferred for the concentrations 
in GMCs {{p/p)m ~ 6 - 8) from ^^CO molecular-line studies (cf. Paper III; Bally et al (1987); 
Williams, Blitz, &; Stark (1995)). The corresponding rms mass-weighted dispersion in p/p\s^ 
3 — 13. Although the density contrast generally increases with the value of the fast-magnetosonic 
Mach number Mp = o"t,/(c^ + (see Fig. 4), there is no obvious one-to-one functional relation 

between M., j3 = (?slv\ q, and the density contrast. In particular, the result obtained by Nordlund 
h Padoan (1999) for purely hydrodynamic quasi-steady turbulence of the relation between the 
density contrast and the Mach number M. does not carry over for (evolving) MHD turbulence. 
When ~ 5 — 9, the Nordlund h Padoan (1999) quasi-steady hydrodynamic-turbulence result 
would predict mass-weighted means and dispersions of pjp in the range ~ 7 — 21 and ~ 18 — 95, 
respectively, larger than the range wc find for our MHD models. Further investigation would be 
required to determine whether, for quasi-steady MHD turbulence, it is possible to find a clean 
functional relation between the mean (and dispersion) oi p/p and the dimensionless parameters 
M. and /? that is independent of the particular instantaneous turbulent power spectrum. Since, 
however, we expect that cold-ISM turbulence is subject to significant transient effects, and in 
addition large "cosmic variance" may result from low-fc dominance of the power spectrum, a one- 
to-one relation of this kind would probably not be realized in GMCs in any case. 
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2. The distributions of column densities N also follow an approximately log-normal form, with 
mean logarithmic contrasts (log(A^/iV)) an order of magnitude smaller that the mean logarithmic 
density contrasts {\og{p/p)) (see Figs. 5-8 and Table 2). The mass-weighted mean column density 
is thus just 7-35% greater than the area- weighted column density N = pL = Mjl? ^ and the mass- 
weighted (area-weighted) dispersion in N is 0.3-0.8 (0.3-0.6) times TV, for models with the range 
of Mach numbers and magnetic field strengths we have analyzed. Large-scale spatial correlation of 
the density perturbations is indicated by the log-normal, rather than Gaussian, form of the column 
density distributions. These large-scale spatial density correlations are associated with large-scale 
correlations in the velocity and magnetic fields (which, for the present models, are input in the 
initial conditions). 

3. We use a binning algorithm to investigate the distribution of kinetic properties for plane-of-sky 
"clumps" as a function of the spatial clump scale. The clumps at a given scale are simply regions 
that contrast with the mean surface density at that scale; Fig. 9 shows an example of the maps 
of these so-defined regions of contrast ("ROCs") at different scales. We find that there is a scatter 
in the values of line-of sight velocity dispersion averaged over projected regions of area s^, ranging 
from the mean value of the velocity dispersion averaged over all cubes of size s'^ up to the velocity 
dispersion for the entire simulation of size 1? (see Figs. 10, 11). The large velocity dispersions 
( "linewidths" ) arise due to the superposition of many small volume elements with differing mean 
velocity along lines of sight through a given projected area. We show that ROCs which have single- 
component velocity distributions ("line profiles") generally consist of several spatially separated 
components along the line-of-sight (see Figs. 12-16). Thus, what looks like a "clump" either on 
the plane of the sky or in position-velocity datacubc may in fact be a superposition of spatially 
unconnected parts with smoothly overlapping velocity distributions. 

4. We find that when mean magnetic fields Bq are weak, there is a large scatter in the distribution 
of total magnetic field strength |B| in any given density regime, and that the mean value of |B| 
varies strongly with density (Fig. 18). For strong mean magnetic fields i?o at moderate densities, 
there is less scatter in the distribution of |B|, and a weaker variation of the mean value of |B| 
with density. At high densities, the variation of the mean of |B| with density is similar for all the 
models (Fig. 19), and is comparable to the indications of increasing field strength from Zeeman 
observations in high density tracers (Crutcher 1999). 

5. We show that, for the models considered here, the line-of-sight average of the line-of-sight 

magnetic field, {Bioi^-, can vary widely across a projected map, and is not positively correlated with 
column density (see Figs. 20, 21). Because of the large scatter in (Bios), an accurate observational 
determination of the global average of {Bios) from the Zeeman method might require a very large 
number of pointings. 

6. We have created simulated polarized-extinction maps by integrating the Stokes parameters along 
lines of sight for different simulation snapshots and orientations (Figs. 22, 23). Because models with 
weak Bq have much more variation in the vector direction of B at a given Mach number, the result 
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is that they have lower average values of the polarization and larger dispersions in the polarization 

angle than their high-i?o counterparts (Fig. 24). We show that the Chandrasekhar- Fermi formula 
can be slightly modified to estimate the plane-of-sky magnetic field strength Bp in terms of the 
mean density p, line-of-sight velocity dispersion 6vios, and plane of sky polarization angle 6(p as 
Bp « 1.8p^/^6vios/d(l), provided 5(j) < 25°. 

One finding particularly notable for the interpretation of molecular-line observations concerns 
the nature of what have long been thought of as large-scale turbulent "clumps" within GMCs. The 
large scatter and relatively shallow mean slope for linewidth vs. projected size distributions of 
observed "clumps" (e.g. Williams, de Geus, &; Blitz (1994); Stutzki &; Glisten (1990)) is similar 
to that found for projected "regions of contrast" (ROCs) in our simulated clouds. The observed 
clumps are identified as coherent structures (peaks and their surroundings) in position-velocity 
molecular line data cubes; our ROCs are overdense on the plane of the sky and generally have 
single-component linewidths. 

Because their linewidths are comparable to those of the larger cloud complexes in which they 
reside, the coherent structures in observed l — b—v cubes have been interpreted as pressure-confined 
clumps (Bertoldi &; McKee (1992)). We believe that these apparent pressure-confined clumps may 
in fact often be line-of-sight superpositions of spatially-unconnected condensations that collectively 
sample from the full turbulent velocity dispersion along the line-of-sight - hence appearing to have 
internal pressure comparable to that of the parent cloud as a whole. We find that the lower envelope 
of our linewidth-size distribution for ROCs closely follows the relation between mean linewidth and 
3D cell size; it may be possible to apply a similar binning procedure to observational maps in order 
to deduce the true 3D mean linewidth-size relation. 

Our structural analysis supports a revision (see references in §1 and §5 for related work) in 
the understanding of GMC clumps that parallels the recent paradigm shift in interpreting the 
cosmological Lyman-alpha forest. Namely, we suggest that the linewidths of these apparent clumps 
is due not to small-scale supersonic "microturbulcncc" , but to a superposition of bulk fiows with 
a large range of correlation lengths. Unlike the situation in cosmology, however, we do not have 
a large-scale Hubble fiow to spread out the velocity field and help us distinguish foreground and 
background concentrations. As a consequence, what appears to be a clump in projection may not, 
in three dimensions, be spatially compact or connected at all! 

A further difference with the cosmological situation is that the overdense regions in clouds are 
not, in general, fated to collapse. The overdense regions are transient objects that form, and then 
disperse, from the effects of time-dependent velocity and magnetic fields. Eventually, through the 
overall dissipation of turbulence and the random superposition of temporary concentrations, a frac- 
tion of the material in a cloud must reach high enough densities to become strongly self-gravitating. 
The subsequent collapse and fragmentation, perhaps initiated at many independent sites in a cloud, 
must ultimately produce a collection of stars. The next generation of numerical simulations will 
require adaptive mesh algorithms to follow this gravitational collapse and fragmentation. Crucial 
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questions are whether the spectrum of stellar masses that forms can ultimately be traced back 
to the bulk initial conditions (e.g. mean density, temperature, magnetic field strength, velocity 
dispersion) in the parent turbulent cloud, and what factors determine the overall rate of conversion 
of gas to stars. 
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Fig. 1. — Energy evolution of model clouds. Models with (5 = 0.01, 0.1, and 1 are shown with 
dotted, solid, and dashed curves, respectively. 
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Fig. 2. — Evolution of turbulent energy E^urb = + ^Eb in model clouds. Models with /3 = 0.01, 
0.1, and 1 are shown with dotted, solid, and dashed curves, respectively. The dot-dash lines indicate 
a slope (xt~^, for reference. 
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Fig. 3. — Comparative statistics of volume density in three model snapshots (B2, C2, D2 from 
Table 2) with matched Mach numbers M. ^ 1 . Solid curves show fraction of volume as a function 
of density relative to the mean (p/p); dashed curves show fraction of mass as a function of p/p. 
Dotted curves show lognormal distributions with the same mean and dispersion as in each model 
snapshot. 
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Fig. 4. — Mass- and volume- weighted mean values of the logarithmic density contrast as a func- 
tion of mean square fast magnetosonic Mach number from model snapshots in Table 2. Pen- 
tagons, squares, and triangles correspond to /3 = 0.01,0.1, and 1, respectively. Dashed lines show 
(log(p/p))M,y = ±0.2(log(l + Ml) + 1) 
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Fig. 5. — Comparative statistics of column density in three model snapshots (B2, C2, D2 from Table 
2) with matched Mach numbers. Projection is along the z axis (perpendicular to the magnetic field). 
In each frame, left-displaced curves show fraction of projected area as a function of column density 
relative to the mean {N/pL = N/N); right-displaced curves show fraction of mass as a function of 
N/N. Dotted curves show lognormal distributions with the same mean and dispersion as in each 
model snapshot. 




Fig. 6. — Comparative statistics of column density at full simulation resolution (a) and at resolution 
a factor four larger in linear scale (b); simulation data is from snapshot B2. 
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Fig. 7. — Same as in 5, except line-of-sight integration is only over z > L/2. 
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Fig. 8. — Comparison of (a) cumulative column density distribution for simulated cloud model (data 
from Fig. 7, with hexagons, squares, and triangles marking the /3 = 1, 0.1, 0.01 model distributions) 
and (b) cumulative extinction distribution for the cloud IC5146. Dotted curves in both panels show 
cumulative log-normal fits. 
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Fig. 9. — Identification of regions of contrast (ROCs) as a function of spatial scale for data from 
model snapshot B2. Panels (a)-(g) outline regions that meet the contrast criterion at increasingly 
fine spatial resolution. A greyscale representation of the projected density is shown in (h) for 
comparison. 
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Fig. 10. — Scale dependence of kinetic quantities in a projected map. In (a)-(d), each point 
represents one of the square regions of contrast (ROCs) identified in Fig. (9), with edge size s. In 
(a), we plot vs. size s (normalized to the box length L) the dispersion of the line-of-sight velocity 
ay in each projected square ROC, and also show (solid line) the mean dispersion in line-of-sight 
velocity for 3D cubes of side ,s. In (b), wc plot vs. s/L the mass M of each ROC (normalized 
to the total mass in the simulation box pL^. In (c), we show vs. M. In (d), we show the 
virial parameter a vs. M. In each frame, dashed lines represent linear least-squares fits to the 
data; dotted lines represent 1-sigma deviations from the fit. In (c) and (d), we plot points from 
different-sized regions with different expansion factors. 
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Fig. 11.— 



Same as in Fig. 10, for model snapshot D2. 
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Fig. 12. — Distribution of mass with line-of-sight velocity for model snapshot B2 projected in the 
z direction, for regions of linear size L/8. Regions meeting surface density contrast criterion (cf. 
Fig. 9) are indicated with heavy outlines. 
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Fig. 13. — Distribution of mass with line-of-sight position for model snapshot B2 projected in the 
z direction, for regions of hnear size L/8. Regions meeting surface density contrast criterion (cf. 
Fig. 9) are indicated with heavy outhnes. 
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Fig. 14. — Same as Fig. (12), for projected region size s = L/16. 
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Fig. 15.— 



Same as Fig. (13), for projected region size s = L/16. 
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Fig. 16. — Same as Figs. (12), (13), for Pmin = lOp- 



line-of-sighl position 
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Fig. 17. — Comparative statistics of magnetic field components for three model snapshots (B2, C2, 
D2 from Table 2) with matched Mach numbers. Solid, dotted, and dashed curves show fraction of 
mass as a function of B^, By, and Bz, respectively. Mean field strengths are nonzero only in the x 
direction, with Bo/(47rpc^)^/^ = 10, 3.16, and 1, for the top, middle, and bottom panels. 
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Fig. 18. — Dependence of the total magnetic field strength on density in three model snapshots 
(B2, C2, D2 from Table 2) with matched Mach numbers. Triangles, squares, and pentagons show 
the mean value of 5^ in each density bin for Mach-7 models with j3 = 0.01, 0.1, and 1, respectively. 
Shaded regions surrounding each curve corresponds to the 1 — cr departures from the mean (errors 
in the means from counting statistics are smaller than the symbols shown). Dashed horizontal lines 
show the values of the square of the mean magnetic field, for the three models. Left and bottom 
scales give magnetic field strength and density p in dimensionless units; right and top scales 
give corresponding fiducial values of and raj^j assuming T = lOK and nj^j = lOOcm"^ for the 
temperature and volume-averaged density. 
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Fig. 19. — Density vs. mean total magnetic field strength at high densities for Mach-9 (top row) 
and Mach-7 (bottom row) model snapshots (fiducial values assume T = lOK and njjj = lOOcm"^). 
Error bars show expected (Poisson noise) error in determination of the means. Also shown (solid 
lines) are fits to \og{B) vs. log(n) for njjj > lO^cm"^; the corresponding slope is indicated in each 
panel. 
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Fig. 20. — Mean line-of-sight-averaged magnetic field strength for Mach-7 model snapshots (B2 — 
left column, C2 - center column, D2 - right column) viewed from three orientations: along the mean 
field (x; top row), at 45° relative to the mean field (center row), perpendicular to the mean field 
{z, bottom row). Contributions to < Bi^g > are weighted by density, and include only zones with 
p/p > 10. Point size is scaled linearly by the value of | < Bigg > |, with positive and negative values 
as shown by the key. For fiducial dimensional B values, we adopt T = lOK and = lOOcm"^. 
Only every 15th point in each direction on the grid is plotted, for clarity. Greyscale underlay shows 
the total column density for each projection. 
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Fig. 21. — Mean line-of-sight-averaged magnetic field strength {Biog) vs. column density N of gas 
at p/p > 10 for model projections shown in Fig. 20. All points at N/pL > 0.5 are plotted; straight 
solid lines show linear least-squares fits. Dashed horizontal lines show the value of the volume- 
averaged mean line-of-sight field Bq sin i for each projection (i is the angle between the plane of sky 
and Bo). 
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Fig. 22. — Column density (color scale, with units pL) and simulated polarization map for model 
snapshot B2 (/3 = 0.01, = 7), projected along z perpendicular to the mean magnetic field. 
The fractional polarization at each point is proportional to the value of a fiducial polarization P 
corresponding to a uniform medium and uniform magnetic field perpendicular to the line of sight, 
arbitrarily set here to P = 0.1 as shown in the key. 




Fig. 23. — Column density and polarization map for model snapshot D2 {(3 = 1, Ai = 7), projected 
along z; definitions as in Fig. 22. 
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Fig. 24. — Distribution of polarization angles for projections of models with matched kinetic energy 
and differing mean magnetic field strength. Left, center, and right columns respectively show 
projections of strong, moderate, and weak-i?o cases (from M = 7 snapshots B2, C2, D2), for the 
mean magnetic field direction lying at varying angles with respect to the plane of the sky (0, 30, 45, 
60, and 90°, from top to bottom). Each panel shows a histogram of the distribution of polarization 
position angles with respect to the most-frequent direction, in degrees. Labels in each panel give 
the dispersion of the distribution. Dashed curves show Gaussian fits. 



